home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_09_08 / 9n08056a < prev    next >
Text File  |  1991-06-19  |  2KB  |  65 lines
  1. /***************** Listing 3 *************************
  2.  
  3.             Cubic Spline Interpolation
  4.  
  5.    double cubic(image, x, y)
  6.      float image[4][4], x, y;
  7.  
  8.    image: pointer to the four values of grid i.e.
  9.           a 4 x 4 array ( image[4][4] )
  10.    x    : sample coordinate
  11.    y    : line coordinate
  12.  
  13.           (0,0)           (0,3)   
  14.             *----*----*----*
  15.             |    |    |    |
  16.             *----*----*----*
  17.             |    | o  |    |
  18.             *----*----*----*
  19.             |    |    |    | 
  20.             *----*----*----*
  21.           (3,0)           (3,3)
  22.  
  23.     If point "o" falls on upper left corner, then 
  24.     return image(1,1). The point should never fall on 
  25.     any of the other corners because the calling 
  26.     program will ensure against this according to the
  27.     filling sequence of the array image[4][4].
  28.  
  29.     Interpolation is first performed in the sample
  30.     direction and then in the line direction, 
  31.     using separability.
  32.  
  33. ******************************************************
  34. #include <stdio.h>
  35.  
  36. double cubic(image, x, y)
  37.  float image[4][4], x, y;
  38. {
  39.   register i;
  40.   float p, q;
  41.   double I_x[4];
  42.   double val = 0.0;
  43.  
  44.   double a1, a2, a3, a4;
  45.   double b1, b2, b3, b4;
  46.  
  47.   p = x - (int) x;
  48.   q = y - (int) y;
  49.   if( (p == 0.0) && (q == 0.0) )
  50.      return( (double) image[1][1] ); 
  51.  
  52.   a1 = p*(1-p)(1-p);      b1 = q*(1-q)(1-q);
  53.   a2 = 1 - 2*p*p + p*p*p; b2 = 1 - 2*q*q + q*q*q;
  54.   a3 = p*( 1 + p - p*p ); b3 = q*( 1 + q - q*q );
  55.   a4 = p*p*(1-p);         b4 = q*q*(1-q);
  56.  
  57.   for(i=0 ; i < 4 ; i++)
  58.     I_x[i] = -a1*image[i][0] + a2*image[i][1] +      \
  59.                 a3*image[i][2] - a4*image[i][3] ;
  60.  
  61.   val = -b1*I_x[0] + b2*I_x[1] + b3*I_x[2] - b4*I_x[3];
  62.  
  63.   return( val );
  64. }
  65.